using System;
using System.Collections.Generic;
using System.Diagnostics;

////////////////////////////////////////////////////////////////////////////////
///
/// Beats-per-minute (BPM) detection routine.
///
/// The beat detection algorithm works as follows:
/// - Use function 'inputSamples' to input a chunks of samples to the class for
///   analysis. It's a good idea to enter a large sound file or stream in smallish
///   chunks of around few kilosamples in order not to extinguish too much RAM memory.
/// - Inputted sound data is decimated to approx 500 Hz to reduce calculation burden,
///   which is basically ok as low (bass) frequencies mostly determine the beat rate.
///   Simple averaging is used for anti-alias filtering because the resulting signal
///   quality isn't of that high importance.
/// - Decimated sound data is enveloped, i.e. the amplitude shape is detected by
///   taking absolute value that's smoothed by sliding average. Signal levels that
///   are below a couple of times the general RMS amplitude level are cut away to
///   leave only notable peaks there.
/// - Repeating sound patterns (e.g. beats) are detected by calculating short-term
///   autocorrelation function of the enveloped signal.
/// - After whole sound data file has been analyzed as above, the bpm level is
///   detected by function 'getBpm' that finds the highest peak of the autocorrelation
///   function, calculates it's precise location and converts this reading to bpm's.
///
/// Author        : Copyright (c) Olli Parviainen
/// Author e-mail : oparviai 'at' iki.fi
/// SoundTouch WWW: http://www.surina.net/soundtouch
///
////////////////////////////////////////////////////////////////////////////////
//
// Last changed  : $Date: 2009-02-21 18:00:14 +0200 (Sat, 21 Feb 2009) $
// File revision : $Revision: 4 $
//
// $Id: BPMDetect.cpp 63 2009-02-21 16:00:14Z oparviai $
//
////////////////////////////////////////////////////////////////////////////////
//
// License :
//
//  SoundTouch audio processing library
//  Copyright (c) Olli Parviainen
//
//  This library is free software; you can redistribute it and/or
//  modify it under the terms of the GNU Lesser General Public
//  License as published by the Free Software Foundation; either
//  version 2.1 of the License, or (at your option) any later version.
//
//  This library is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//  Lesser General Public License for more details.
//
//  You should have received a copy of the GNU Lesser General Public
//  License along with this library; if not, write to the Free Software
//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
//
////////////////////////////////////////////////////////////////////////////////

namespace soundtouch
{
	public class BPMDetect
	{
		const int INPUT_BLOCK_SAMPLES = 2048;
		const int DECIMATED_BLOCK_SAMPLES = 256;
		const int MIN_BPM = 29;
		const int MAX_BPM = 230;
		
		/// Auto-correlation accumulator bins.
		float[] xcorr;

		/// Amplitude envelope sliding average approximation level accumulator
		float envelopeAccu;

		/// RMS volume sliding average approximation level accumulator
		float RMSVolumeAccu;

		/// Sample average counter.
		int decimateCount;

		/// Sample average accumulator for FIFO-like decimation.
		float decimateSum;

		/// Decimate sound by this coefficient to reach approx. 500 Hz.
		int decimateBy;

		/// Auto-correlation window length
		int windowLen;

		/// Number of channels (1 = mono, 2 = stereo)
		int channels;

		/// sample rate
		int sampleRate;

		/// Beginning of auto-correlation window: Autocorrelation isn't being updated for
		/// the first these many correlation bins.
		int windowStart;

		/// decay constant for calculating RMS volume sliding average approximation
		/// (time constant is about 10 sec)
		const float avgdecay = 0.99986f;

		/// Normalization coefficient for calculating RMS sliding average approximation.
		const float avgnorm = (1 - avgdecay);
		
		/// FIFO-buffer for decimated processing samples.
    	FIFOSampleBuffer buffer;

		public BPMDetect(int numChannels, int aSampleRate)
		{
			this.sampleRate = aSampleRate;
			this.channels = numChannels;

			decimateSum = 0;
			decimateCount = 0;

			envelopeAccu = 0;

			// Initialize RMS volume accumulator to RMS level of 3000 (out of 32768) that's
			// a typical RMS signal level value for song data. This value is then adapted
			// to the actual level during processing.
			RMSVolumeAccu = (0.092f * 0.092f) / avgnorm;

			// choose decimation factor so that result is approx. 500 Hz
			decimateBy = sampleRate / 500;
			Debug.Assert(decimateBy > 0);
			Debug.Assert(INPUT_BLOCK_SAMPLES < decimateBy * DECIMATED_BLOCK_SAMPLES);

			// Calculate window length & starting item according to desired min & max bpms
			windowLen = (60 * sampleRate) / (decimateBy * MIN_BPM);
			windowStart = (60 * sampleRate) / (decimateBy * MAX_BPM);

			Debug.Assert(windowLen > windowStart);

			// allocate new working objects
			xcorr = new float[windowLen];

			// allocate processing buffer
			buffer = new FIFOSampleBuffer();
			// we do processing in mono mode
			buffer.setChannels(1);
			buffer.clear();
		}

		/// convert to mono, low-pass filter & decimate to about 500 Hz.
		/// return number of outputted samples.
		///
		/// Decimation is used to remove the unnecessary frequencies and thus to reduce
		/// the amount of data needed to be processed as calculating autocorrelation
		/// function is a very-very heavy operation.
		///
		/// Anti-alias filtering is done simply by averaging the samples. This is really a
		/// poor-man's anti-alias filtering, but it's not so critical in this kind of application
		/// (it'd also be difficult to design a high-quality filter with steep cut-off at very
		/// narrow band)
		private uint decimate(float[] dest, float[] src, uint  startIndex, uint numsamples)
		{
			uint count, outcount;
			double output;

			Debug.Assert(channels > 0);
			Debug.Assert(decimateBy > 0);
			outcount = 0;
			for (count = startIndex; count < startIndex + numsamples; count ++)
			{
				uint j;

				// convert to mono and accumulate
				for (j = count; j < count + channels; j ++)
				{
					decimateSum += src[j];
				}

				decimateCount ++;
				if (decimateCount >= decimateBy)
				{
					// Store every Nth sample only
					output = (double)(decimateSum / (decimateBy * channels));
					decimateSum = 0;
					decimateCount = 0;
					dest[outcount] = (float)output;
					outcount ++;
				}
			}
			return outcount;
		}



		// Calculates autocorrelation function of the sample history buffer
		void updateXCorr(uint process_samples)
		{
			int offs;
			List<float> pBuffer;

			Debug.Assert(buffer.numSamples() >= (uint)(process_samples + windowLen));

			pBuffer = this.buffer.buffer;
			for (offs = windowStart; offs < windowLen; offs ++)
			{
				double sum;
				int i;

				sum = 0;
				for (i = 0; i < process_samples; i ++)
				{
					sum += pBuffer[i] * pBuffer[i + offs];    // scaling the sub-result shouldn't be necessary
				}
				//        xcorr[offs] *= xcorr_decay;   // decay 'xcorr' here with suitable coefficients
				// if it's desired that the system adapts automatically to
				// various bpms, e.g. in processing continouos music stream.
				// The 'xcorr_decay' should be a value that's smaller than but
				// close to one, and should also depend on 'process_samples' value.

				xcorr[offs] += (float)sum;
			}
		}



		// Calculates envelope of the sample data
		private void calcEnvelope(float[] samples, uint numsamples)
		{
			const float decay = 0.7f;               // decay constant for smoothing the envelope
			const float norm = (1 - decay);

			uint i;
			double output;
			float val;

			for (i = 0; i < numsamples; i ++)
			{
				// calc average RMS volume
				RMSVolumeAccu *= avgdecay;
				val = (float)Math.Abs((float)samples[i]);
				RMSVolumeAccu += val * val;

				// cut amplitudes that are below 2 times average RMS volume
				// (we're interested in peak values, not the silent moments)
				val -= 2 * (float)Math.Sqrt(RMSVolumeAccu * avgnorm);
				val = (val > 0) ? val : 0;

				// smooth amplitude envelope
				envelopeAccu *= decay;
				envelopeAccu += val;
				output = (double)(envelopeAccu * norm);

				samples[i] = (float)output;
			}
		}

		public void inputSamples(float[] samples, uint numSamples)
		{
			float[] decimated = new float[DECIMATED_BLOCK_SAMPLES];

			// iterate so that max INPUT_BLOCK_SAMPLES processed per iteration
			while (numSamples > 0)
			{
				uint block;
				uint decSamples;
				uint startIndex = 0;

				block = (numSamples > INPUT_BLOCK_SAMPLES) ? INPUT_BLOCK_SAMPLES : numSamples;

				// decimate. note that converts to mono at the same time
				decSamples = decimate(decimated, samples, startIndex, block);
				startIndex += Convert.ToUInt32(block * channels);
				numSamples -= block;

				// envelope new samples and add them to buffer
				calcEnvelope(decimated, decSamples);
				buffer.putSamples(decimated, decSamples);
			}

			// when the buffer has enought samples for processing...
			if ((int)buffer.numSamples() > windowLen)
			{
				uint processLength;

				// how many samples are processed
				processLength = Convert.ToUInt32(buffer.numSamples() - windowLen);

				// ... calculate autocorrelations for oldest samples...
				updateXCorr(processLength);
				// ... and remove them from the buffer
				buffer.receiveSamples(processLength);
			}
		}

		public float getBpm()
		{
			double peakPos;
			PeakFinder peakFinder = new PeakFinder();

			// find peak position
			peakPos = peakFinder.detectPeak(xcorr, windowStart, windowLen);

			Debug.Assert(decimateBy != 0);
			if (peakPos < 1e-6) return 0.0f; // detection failed.

			// calculate BPM
			return (float)(60.0 * (((double)sampleRate / (double)decimateBy) / peakPos));
		}
		
	}
}


